function Z = pp_baugh_wooley(A, B)
    len = length(A);
    Z = zeros(len+1, length(A)+len);
    C = zeros(len+1, length(A)+len);
    for i = 1:len,
        P = A(len-(i-1))*B;
        s = P(1);
        
        C(i,:) = [zeros(1,len*2-len+1-i) 1 zeros(1,len-2+i)];
        
        P(1) = not(P(1));
        if(i == len)
            P = not(P);
        end
        P = [ones(1,len-i+1)*s P zeros(1,i-1)];
        Z(i,:) = P;
    end
    C = mult_helper(C(1:len,:));
    C = [not(C(1:len+1)) zeros(1,len-1)];
    D = [zeros(1,len) 1 zeros(1,len-1)];
    C = add_sub(C,D,0);
    E = [zeros(1,i) 1 zeros(1,i-1)];
    C = add_sub(C,E,0);
    Z(i+1,:) = C;
end

function [S ci] = add_sub(A,B,ci)
    len = length(A);
    sub = ci;
    for i = len:-1:1,
        [S(i),ci] = bit_add(A(i),xor(B(i),sub),ci);
    end
end

function [s,co] = bit_add(a,b,ci)
    t = xor(a,b);
    s = xor(t,ci); 
    co = (a&&b)||(t&&ci);
end